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ABSTRACT 



Context. Long-term precise Doppler measurements with the CORALIE spectrograph revealed the presence of two massive compan- 
ions to the solar-type star HD202206. Although the three-body fit of the system is unstable, it was shown that a 5:1 mean motion 
■ resonance exists close to the best fit, where the system is stable. It was also hinted that stable solutions with a wide range of mutual 

inclinations and low 0-C were possible. 

Aims. We present here an extensive dynamical study of the HD202206 system aiming at constraining the inclinations of the two 
(~ i , known companions, from which we derive possible ranges of value for the companion masses. 

^ |. Methods. We consider each inclination and one of the longitude of ascending node as free parameters. For any chosen triplet of these 

I ' parameters, we compute a new fit. Then we study the long term stability in a small (in terms of O-C) neighborhood using Laskar's 

O ' frequency map analysis. We also introduce a numerical method based on frequency analysis to determine the center of libration mode 

inside a mean motion resonance. 

^ . Results. We find that acceptable coplanar configurations (with low stable orbits) are limited to inclinations to the line of sight 

^ ' between 30° and 90°. This limits the masses of t50th companions to roughly twice the minimum: 6 [16.6 Mj„,,; 33.5 Mj„,,] and 

6 [2.2Mj„p; 4.4Mj„p]. Non coplanar configurations are possible for a wide range of mutual inclinations from 0° to 90°, although 
AO = 0[;r] configurations seem to be favored. We also confirm the 5:1 mean motion resonance to be most likely. In the coplanar 
^ ■ edge-on case, we provide a very good stable solution in the resonance, whose x~ does not differ significantly from the best fit. Using 

• our method to determine the center of libration, we further refine this solution to obtain an orbit with a very low amplitude of libration, 

' as we expect dissipative effects to have dampened the libration. 

Key words. HD202206 - extra-solar planetary systems - mean motion resonance - frequency map analysis 
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1. Introduction companion made the second planet evolution very chaotic, ex- 

cept for initial conditions in the 5:1 mean motion resonance. 



O^^v^jjL iwi 1111 LI ai \^wii\aiLi\Jiia 111 Liiv^ ^ . 1 niveau iiiwLiwii iv^awiiaiiv^v^. 

The CORALIE planet-search program in the southern hemi- ^-^^^ ^j^^ associated resonant island actually lies close to the 

^ sphere has found two companions around the HD202206 star j^inimum value of the best fit, they concluded that the system 

• rH . The first one is a very massive body with 17.5 Mj„p imm- ^j^^^j^ j^^j^^^ i„ ^^^^ 5.^ resonance. Later on, Gozdziewski 

p<; ■ mum mass (Vdry et al."2002'), while the second companion is co-workers also looked for stable solutions in this system us- 

a 2.4 Mjup minimum mass planet dCorreia et al. 2005). The par- -^^ ^^^^-^ Q^^p algorithm (.Gozdziewski et al.ll2006) . They pro- 

ent star has a mass of 1 .044 solar masses, and is located 46.3 pc ^j^^^ p^^^i^le solutions (among many others), one coplanar, 

from the Solar System. The HD202206 planetary system is an ^^^^ ^ ^^^y high mutual inclination. 

interesting case to investigate the brown dwarf desert since the „. ^, j . u 1. • j • 

°. , , ,j. , . Since then, new data have been acquired using the 

more massive companion can be either a huge planet (formed in ^^^r. a ? ti- . u a j j . u j 

, . , , ■ ^ ,• I CORALIE spectrograph. A new reduction of the data changed 

the circumstellar disk) or a low-mass brown dwarf candidate. ^ . . .. ^, „ ^, m^^n^^n^ , » 

r- : 71 rrTTT-^i ^ , , , , • , , some parameters, including the mass of the HD202206 star A 

Correia et alj (K005D round that the orbital parameters ob- r-. ■ ^ a c a ■ a 

. ; — — ^ — - — , . . . ^ , . new fit assuming a coplanar edge-on configuration was derived 

tamed with best fit for the two planets leads to a catastrophic ^^^^ ^^^.^j ^^j^^j ^^^^ ^^^^ ^^^^^^^^ ^^jH 

events in a short time (two keplerians fit and full three-body fit ^ u • »u c 1 u » • 1 » n 

^ . , , , • , , appears to be in the 5:1 mean motion resonance, but IS also still 

alike). This was not completely unexpected given the very large ui -ru » » 1 • a- a e '■ — mJinnck 

' ... ,„ r r , , unstable. The most striking dinerence trom lCorreia et al.l (I2005D 

eccentricities (0.435 and 0.2b/) and masses of the two plan- ;„ .1,= „f uT-nmon<;„ 

TT • , ■ I t . — N ■ IS thc smallcr eccentiicity ot HD20220o c . 

ets. Using frequency analysis (lLaskaiill990l 1 19931 119991) they ^ ^ , .„ . . ^ •, u j 

f J / J 1 u 1 J ■ J .1. u » £^ In the present work, we will continue in more detail the dy- 

performed a study on the global dynamics around the best fit, ■ , , ■ • it; ■■ n j-o^^ t. • 1 i 

AC A.u ..u . v 1 • . .u c . namical study started in ICorreia et all (I2005D . using the new fit 

and found that the strong gravitational interactions with the first . ■' . ... ' : — , . , ,. , 

as a starting point. We also aim to find constraints on the orbital 

* Based on observations made with the CORALIE instrument on the Parameters of the two known bodies of this system, in partic- 

EULER 1.2m telescope at La Silla Observatory under the?? programme "^ar the inchnations (and thus the real masses of the planets). 

ID ??. The table with the radial velocities is available in electronic form Gozdziewski et al. (2006) already showed that stable fits could 

at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or be obtained with different inchnations, using a particular fitting 

via,http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A-HA/????) genetic algorithm that adds stability computation to select its 
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populations (GAMP). Although very effective in finding a stable 
fit, this algorithm cannot find all possible solutions. We prefer 
here an approach which separates the fitting procedure from dy- 
namical considerations, as it allows for a better assessment of the 
goodness of the fit, and whether the model is a good description 
of the available data. The trade-off is a more difficult handling of 
the high number of parameters. 

We briefly present the new set of data in Sect. |2l the numer- 
ical methodology in Sect. [3] We review in details the dynamics 
in the coplanar edge-on case in Sect.|4] We then release the con- 
straint on the inclination of the system from the line of sight in 
Sect.|5] and finally briefly investigate the mutually inclined con- 
figurations in Sect.m 

2. New orbital solution for the HD202206 system 



Table 1. Best Newtonian fit SI for the HD202206 system as- 
suming ib = i, = 90°, and AQ = 0°. 



Param. 


SI 


inner 


outer 


Vo 


[km/s] 


14.729 


±0.001 


P 


[days] 


256.389 ± 0.044 


1397.445 ± 19.056 


K 


[m/s] 


564.82 ± 1.42 


38.08 ± 1.21 


e 




0.431 ±0.001 


0.104 ±0.024 


O) 


[deg] 


161.91+0.27 


105.56 ± 15.77 


M + cij 


[deg] 


239.016 ±0.13 


250.38 ±2.71 


a 


[AU] 


0.8053 


2.4832 


i 


[deg] 


90 


90 


a 


[deg] 








m 


[Mjup] 


16.59 


2.179 


Date 


[JD-2400000] 


53000.00 


rms 


[m/s] 


7.4544 






1.411 



Errors are given by the standard deviation cr. This fit corresponds to a 
coplanar system seen edge-on, with minimum values for the masses. 



The CORALIE observations of HD202206 started in August 
1999 and the last point acquired in our sample dates from 
September 2006, corresponding to about seven years of obser- 
vations and 92 radial-velocity meas urements. Using the itera- 
tive Levenberg-Marquardt method ( Press et al.lll992h we fit the 
observational data using a 3-body Newtonian model, assuming 
co-planar motion perpendicular to the plane of the sk y, simi- 
larly to what has been done in (Correia et al. 2005; Correi a^t al.l 
l2009bl) . Notice tha t we changed t he ref erence date with respect 
to the solution in ICorreia et al.l (|2005|). The mass of th e star 
has also been updated to 1 .044 jSousa et al.1 l2008l) . This 
fit yields two planets with an adjustment of = 1.41 and 
rms = 7.45 ms slightly above the photon noise of the in- 
strument which is around 6.69 ms~'. We confirm th e already 
detected planets (Udry et al.l l2002HCorreia et al.ll200 ?) with im- 
proved orbital parameters, one at P = 254.8 day, e = 0.431, and 
a minimum mass of 16.6 Mjup, and the other at f = 1397 day, 
e - 0.104, and a minimum mass of 2.18 Mjup (Table[T]). In fig- 
ure [U we plot the observational data superimposed on the best 
fitted solution. 

We also fitted the data with a 3-body Newtonian model for 
which the inclination of the orbital planes, as well as the node of 
the outer planet orbit, were free to vary. We were able to find a 
wide variety of configurations, some with low inclination values 
for one or both planets, that slightly improved our fit to a min- 
imum - 1.30 and rm5 = 7.08 m/s. However, all of these 




51500 52000 52500 53000 53500 54000 

JD-2400000 [days] 

Fig, 1, CORALIE radial velocities for HD 202206 superimposed 
on a 3-body Newtonian orbital solution (Table[T]). 



determinations remain uncertain, and since we also increase the 
number of free parameters by three, we cannot say that there has 
been an improvement with respect to the solution presented in 
Tabled 



3. Numerical set-up 

3.1. Conventions 

In this paper, the subscripts b and c will respectively refer to the 
body with shortest and longest orbital period (inner and outer). 
The initial conditions ill-constrained by the radial velocity data 
are it, ic,^b,^c (the inclinations, and longitude of ascending 
nodes). We are using the observers convention which sets the 
plane of sky as the reference plane (see figHJi. As a consequence, 
the nodal line is in the plane of sky, and has no cinematic impact 
on the radial velocities. From the dynamical point of view, only 
the difference between the two lines of nodes AQ matters. In par- 
ticular the mutual inclination depends on this quantity, through: 



cos / = cos ib cos ie + sin sin i^ cos AQ. 



(1) 



Q/, can thus always be set to 0° in the initial conditions, 
which leads to AQ = Q^., and only three parameters are left free 
(ib, ic, Qc)- They are connected to the three interesting unknowns 
of the system, namely the mutual inclination / (Eq.[Tli, and the 
two planetary masses nib, as: 



1/3 



mbinit + nib) _ 
nit + nib + nic {InGY^^ sin ib 

nit + nib + nic (InCy^^ sin ic 



(2) 



where m, is the star mass, Kcr the amplitude of the radial ve- 
locity variations, f „- the orbital period, e^r the eccentricity, and G 
the gravitational constant. Basically, choosing the values of the 
inclinations is akin to setting the two companions masses. And 
for given values of (ib, ic), we can control the mutual inclination 
with Qc through Eq. [1] 

Note that the denominator on the left hand sides in Eq. |2]is 
the total mass of the system. This term comes from the transfor- 
mation to the barycentric coordinates system. As a consequence, 
the two equations are coupled. Of course, we can usually neglect 
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Fig. 2. Angles defining the orbit's orientation in space. We fol- 
low the observers convention which sets the plane of sky as the 
reference frame, for which the edge-on coplanar configuration 
is ib = ic - 90° and AQ - 0°. The observer is hmited to the 
velocity projected on the k axis. 



Finally we notice that one can restrict one of the two incli- 
nations to the line of sight to ]0; n/l], which will be the case for 
ib- 

3.2. Fitting procedure 

The influence of ib, ic, and Q^. on the radial velocity data is usu- 
ally very small, and the perturbations depending on them have 
very long time scales. This makes any attempt to fit those pa- 
rameters virtually impossible at present. Only with very strong 
mea n motion resonances, such as observed in the GJ 876 sys- 
tem (iLaughlin & Chambersir200lHCorreia et al.ll2009bl) . can one 
hope to fit the inclinations. In this case the mean motion res- 
onance introduces important short time scale terms (compared 
with the precision and time span of the observations). For the 
HD202206 system, the set of radial velocity data does not cover 
a long enough period of time . 

Since the three parameters ib,ic and Q.^ are very poorly con- 
strained by the radial velocity data, we cannot fit the data with 
a model which includes them. Instead, for any chosen {ib,ic,^c) 
set, we co mpute a new best fi t using a Levenberg-Marquard min- 
imization dPress et al.lll992l) and a three-body model, but with 
eleven free parameters: the center of mass velocity y, and for 
each planet, the semi-amplitude of radial velocity K, the period 
P, eccentricity e, mean anomaly M and periastron oj, all given at 
the initial epoch. Throughout the paper the initial conditions are 
given at the same initial epoch To -2 453 000 JD. 

At this point we have a complete description of the system 
with the mass of the hosting star and 14 parameters (7 for each 
planets) as follows: 



the companion masses in this term, and decouple the equations. 
However when we change the inclinations, the planetary masses 
will grow to a point where this approximation is no longer valid. 
Here we always solve the complete equations, regardless of the 
inclinations. 

As long as the companions masses are small compared to the 
primary, they are scaled to 1/ sin; to a good approximation. We 
can define two factors kb - 1/ sin ib and k^ - 1/ sin i^. With mf^ 

and mf^ the minimum masses obtained for the edge-on coplanar 
case (see section|4]and tablelTJ, we can write: 



nib ~ kbin 



(0) 



1 



* sin/i * 



(3) 



■,kcmf^^ 



sm Ir 



►,(0) 



For a given factor k, two values of the inclination are possible: x 
andn — X (where x e [0; 7r/2]). For instance ; = 30° and / = 150° 
give k = 2. 

Additionally, for a given pair {ib,ic), the accessible mutual 
inclinations / are limited. Since the inclinations and ic are an- 
gles between and n excluded, prograde coplanar configurations 
are only possible for ib - ic and Af2 = (Eq.[Tli. Similarly, retro- 
grade coplanar configurations are only possible for ib+ic - and 
AQ - n. This means that kb - kc in both cases. More generally, 
/ > \ib - ic\, and: 



ib + ic ^ ^ 
ib + ic > n 



■ I <ib + ic , 
• I < 2n — ib — ic 



(4) 



For any value of /, two values of AQ. are possible: x and -x, 
since AQ. appears through its cosine in Eq. [T] The extrema are 
obtained for AQ = (maximum), and AQ = n (minimum) 



- 4 chosen : /o- and Qq- 

- and 10 fitted : Kg-, Pa, ^a-, Mg-, and Wo-- 

However, masses and elliptic coordinates are easier to ma- 
nipulate for a dynamical study. Using Eq. |2] we obtain a system 
of two equations with two unknowns (nib and nic) which is easily 
solved with a Newton algorithm. The semi-major axis are then 
obtained using Kepler's third law. 



3.3. Numerical integrations 

For the numerical integrations we use the Newtonian equations 
with secular corrections for the relativity. The Newtonian part 
of the integration is carried_out by the symplectic integrator 
SABAC4 of lLaskar&Robutell fcOOlft with a step size of 0.02 
year. The secular corrections for the r elativity are computed 
from t he perturbation formulae given in iLestrade & BretagnonI 
(Il982h : 



dfl 
df 
dA 
dt 



dm 
"d7 
de 
di 



2na^e sin v 



10^ 
r 



na I an" r 
■ n + ^l-r■ h— + 6- - 20 



-10 



2 2 

a rj- 



ari a 



nai] 



-10 



37 



narj sin v 



narj sin v 



a 

10- -7 

r 



^ + 18^-7 



(5) 
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where ji - Gmt/c^, n is the mean motion, and 77 - Vl - e^. 
These equations are averaged to obtain the following first order 
secular perturbations: 



< a : 

< e . 



<M>-- 



jin 



15 



<m >- 



vr 



fl(l 



(6) 



These corrections are computed every 100 steps, that is every 
two years, with the current values at the given step for e, a and 
n, for each planet. These approximated equations have be en suc- 
cessfully tested by comparison with INPOP jFienga et alEOOSi) . 



3.4. Stability tliresliold 

In order to study the sta bility of a given orbi t, we use Laskar's 
frequency map analysis (lLaskaiill990l Il993h . Using a numeri- 
cal integration of the orbit over a time interval of length T, we 
compute a refined determination (in ° /yr) of the mean motions 
«!, «j obtained over two consecutive time intervals of length 
Ti - T 12. The stability index Di - \n'^ - n\ \ provides a mea- 
sure of the chaotic diffusion of the trajectory. Small values close 
to zero correspond to a regular solution while high values are 
synonymous of strong chaotic motion (Laskar 1993). 

In this paper we look at many different orbits for many dif- 
ferent initial conditions to detect stable regions. This calls for a 
way to automatically calibrate a threshold for stability Dw^. To 
that end we use a second stability index D2. Using the same nu- 
merical integration, we compute two new determinations of the 
mean motion «2 and n'^ over two consecutive time intervals of 
length T2 - Ti/k, where k > 1. In the case of quasi-periodic 
motion, the diffusion should be close to zero but it is limited 
by the precision of the determination of the frequencies. Since 
Di is computed over longer time intervals, the frequencies are 
better determined, and thus, Di should be, on average, smaller 
than D2. On the contrary, for chaotic trajectories, the diffusion 
will increase on average for longer time intervals. We can then 
determine an approximated value D\i^ for which: 



Di <s; Di 



Di>D2, 



(7) 



We will then be assured that the first kind of orbits is in general 
stable, while the latter is considered chaotic. 

This approach is best used statistically over a grid of ini- 
tial conditions, especially when we try to use a small integration 
time. In order to determine Diin, for a particular diffusion grid, 
we look at the distribution of Di < D2 trajectories as a function 
of Di. We actually work with smoothed values Dj and of 
Di and D2 to reduce the influence of the chaotic orbits whose 
mean motion diffusions are small by mere chance. They appear 
as low diffusion orbits inside high diffusion region. The smooth- 
ing function is a simple geometric mean over the closest neigh- 
bors. Other functions, such as a convolution with 2D Gaussian, 
have been tested, but do not yield significantly better results. 

We bin the logDj data in 0.5 wide bins, and compute for 
each bin the percentage of Dj < orbits. Fig. [3] shows a typ- 
ical distribution obtained from a diffusion grid (in this case the 
top panel of figure|6]l. It reproduces the behavior expected from 



Eq. I2] low diffusion orbits tend, in great majority, to have their 
diffusion index diminish when time increases. 



Q 
V 



0.6 



0.4 



0.2 




-6 -5 -4 -3-2-1012 
log D, 

Fig. 3. Distribution of Di < D2 trajectories from the top panel of 
Fig.|6l Each trajectory integrated is binned with respect to its dif- 
fusion index log Di, after the diffusion grid has been smoothed. 
For each bin we compute the proportion of trajectories with a 
decreasing diffusion index over time: Di < D2. The histogram 
shows the results for 0.5 wide bins. As expected from equation|7] 
orbits with a high diffusion index Di, and Di < D2 are nearly in- 
existant (logDi > -1), while we observe the opposite situation 
for low diffusion index orbits (logDi < -3). 



We choose to define Du^ as the Di value for which 99% of 
the trajectories exhibit Di < D2. Graphically, it is the abscissa 
for which the curve in Fig. [3] crosses y = 0.99. In this example 
we get log Aim ~ -2.87. 

The 99% threshold is actually a compromise that works in 
the majority of encountered cases: it minimizes the number of 
orbits wrongly flagged as stable (false stable) or unstable (false 
unstable). This is also approximately the value for which the 
number of false stables and false unstables is equivalent. For per- 
centages lower than 95% the number of false unstables is nearly 
null. It is easy to understand since a higher percentage threshold 
implies a lower value for Diim, which in turns leads to very few 
actually unstable orbits flagged, while we might miss several sta- 
ble ones, and vice-versa. To estimate the number of false stables 
and unstables, we recomputed the diffusion grid on a longer in- 
tegration time, 2 x 40000 years, and use this grid as a reference 
(Fig.©. 

4. Review of the coplanar edge-on case 

4.1. Global dynamics 



Following lCorreia et al.l (l2005h . we will study in more details the 
dynamics in the neighborhood of the 3 -body fit obtained in the 
case of coplanar orbits with sin = sin = 1 (that is, the sys- 
tem seen edge-on). The best fit to the radial velocity data for this 
particular configuration is given in table [T] It is different from 
the solution S4 in iCorreia et alj (I2005h (Table 4) as explained 
in Sect. |2] For the dynamical aspect of the system, the impor- 
tant change is the decrease in planet c's eccentricity. As a con- 
sequence regions outside resonances are expected to be more 
stable, and the environment of the fit should be less chaotic. 
However this new solution is still unstable: the outer planet is 
lost shortly after about 150 millions years. 
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false negatives 







0.5 0.6 0.7 0.8 0.9 1 

tfirestiold 

Fig. 4. Percentage of orbits wrongly flagged stable (false stable) 
or unstable (false unstable). Once is choosen using a given 
percentage threshold (see figure[3]), we have a criterion for stable 
{Di < Diim) and unstable orbits (Di > Du^). We compare the 
results with a diffusion grid computed over a longer time interval 
(in this case 2 x 40000 years) taken as a reference. A few orbits 
deemed stable from the reference diffusion grid were thought to 
be unstable and vice versa. The solid line traces the number of 
those faulty orbits for different values of the threshold. It appears 
to be minimum between 0.95 and 1 . The dotted line traces the 
false positives, and the dashed line, the false negatives (the solid 
one is the sum of those two). 



We look for possible nearby stable zones, keeping 
HD202206b parameters constant since they are much better 
constrained, with small standard eiTors. We assume for now that 
the system is coplanar and seen edge-on, that is with both in- 
clination at 90°, and - 0°. We let Oc, ec, and Wc vary. We 
always keep Mc + <±>c constant as it is much better constrained 
by the radial velocity data. This implies that when we change 
lOc, the mean anomaly Mc varies accordingly. In the particular 
case where = 0°, this means that the initial mean longitude 
Ac is kept constant. For each initial conditions we compute the 
diffusion index log D\, and the square root of the reduced 

Fig. |5] shows a global picture of the dynamics around the 
fit, in the planes {qc, e^) and {cic, tOc) of initial conditions. The 
step-sizes for Of, w^, are respectively 0.005AC/, 2° and 0.004. 
The other parameters were kept constant and taken from the fit 

SI (Table[T]i. The level curves give the ^[)^ value computed for 
each set of initial conditions. The color scale gives the diffusion 
index logDi. The yellow to red areas are very chaotic, mainly 
due to the large eccentricities and masses of both planets. 

The orbital solution SI lies inside the ^[)^ - 1.5 level 
curves, at the coordinates marked by a cross, inside a high dif- 
fusion (green) area. Several low diffusion (blue and dark blue) 
zones exist for which the orbits are stabilized either by mean 
motion resonances, or by locking of Am around 0°. Orbits stabi- 
lized by the corotation of the apsidal lines are the blue to black 
zones around Wc - 190° (bottom panel). The width (in the coc di- 
rection) increases with the semi-major axis Qc, from 90 degrees 
at 2 AU, to nearly 360 degrees at 4 AU, since a wider libration of 
Act around 0° is possible without close encounters when the dis- 
tance between the two planets increases. The red to green more 
or less vertical stripes cutting through these zones mark mean 
motion resonances, which for the most part have a destabiliz- 




2.0 2.5 3.0 3.5 4.0 




-6 -4-2 2 

log(D,) 



Fig. 5. Global view of the dynamics of HD202206 for variation 
of the semi-major axis and periastrum of the outer planet (bot- 
tom panel) or semi-major axis and eccentricity (top panel). The 
step-sizes for Gc, (±>c, are respectively Q.QQ5AU, 2° and 0.004. 
The other parameters were kept constant and taken from the fit 
SI (Table [T]i. The color scale is the stability index logDi ob- 
tained through a frequency analysis of the longitude of the outer 
planet over two consecutive time intervals of 8000 years. The 

level curves give the value computed for each choice of pa- 
rameters. The two horizontal black lines mark the intersection of 
the two grids. Most of the orbits are chaotic (yellow to red dots), 
but regular orbits zones exist (blue dots). One of these regular or- 
bits regions lies inside the low yjx^ region: the dark area around 
Qc - 2.5 AU and Uc - 50°. It corresponds to the stable island of 
the 1/5 mean motion resonance. The cross marks the SI orbital 
solution. 



ing effect in the two plans considered in Fig. |5] However, the 
stronger ones also have stable orbits: 

- the 1 /4 MMR at 2.2 AU with two stable islands around (jJc - 
260° andw, = 50°. 

- the 5:1 MMR at 2.5 AU with a notable stable island around 
LOc - 70° where the best fit is located. 

- the 1 /6 MMR at 2.8 AU with a stable island around cDc - 0° 
or high eccentricity. 
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4.2. Stable fit 

We now take a closer look at the 5:1 mean motion resonance 
island around fl^. = 2.5 AU and w^, = 50°, where we believe the 
system is presently locked. Fig.|6]was constructed the same way 
as Fig.|5] but the step-sizes are now 0.0025 AU for Oc, 0.5° for 
ojc, and 0.002 for e^- 

For eccentricities higher than 0.2, orbits are very chaotic (red 
dots), as the outer planet undergoes close encounters with the in- 
ner body. At lower eccentricities we notice some lower diffusion 
orbits for cjc > 120°. Those orbits lie far outside the resonance, 
but may be stable because of the low eccentricity of the outer 
planet, and apsidal locking mechanism. They are however too 

far from the best fit ( > 1 .5) and are less likely to be a good 
guess of the actual configuration of HD202206 system. 

A very noticeable feature of this resonant island appearing in 
both panel, is the existence of two distinct stable regions, sepa- 
rated by chaotic orbits inside the resonance itself. The two stable 
regions actually correspond to two different critical arguments: 
Ab-SAc+A-vJc in the structure on the rim, and Ab-5Ac + mh + 3mc 
in the center. 

The orbital solution SI (white cross) lies very close to the 
latter, in a chaotic region (green and yellow dots) between the 
two stable parts of the resonant island. In fact, one can pick sta- 
ble orbits with a ^/x^ not significantly worse than the best fit. For 
instance, the orbital solution S2 given in table|2]is stable and has 

a of 1.4136 (marked by a filled white circle). The orbital 
elements are the same as SI, except for a^. which was adjusted 
from 2.4832 AU to 2.49 AU. 

Table 2. Stable orbital parameters S2 for the HD202206 system 
for //, = ic = 90° and AO. = 0°. Using the fit SI as a starting 
point (Table [1]), we select a value for the semi-major axis of the 
outer planet such that the system becomes stabilized in the 5:1 
mean motion resonance. This orbit is marked by a white filled 
circle in Fig.|6] 
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4.3. Resonant and secular dynamics 

The orbital solution S2 was integrated over 5 Gyr. It remained 
stable, and displays a regular behavior during the whole time. 

Using frequency analysis on an integration over IMyr, we 
determined its fundamental frequencies (tableO. Following our 
notation, rit, and ric are the mean motions. The secular frequen- 
cies are noted gi and g2. Finally Ig is the frequency associated 
with the resonance's critical angle 6. The fundamental secular 
frequencies gi and g2, related to the periastron of the inner and 
outer planet, correspond to the periods Pi » 14.1 x 10^ yr and 
P2 a: -791 yr (the periastron of the outer planet is retrograde). 

Due to the mean motion resonance, a linear relation links the 
first four fundamental frequencies in Table [3] «/, - Sric + g\ + 




2.40 2.45 2.50 2.55 2.60 




2.40 2.45 2.50 2.55 2.60 



-6 -4-2 2 

log(D,) 

Fig. 6. Global view of the dynamics of HD202206 for variation 
of the semi-major axis and periastrum (bottom), and semi-major 
axis and eccentricity (top) of the outer planet. The step-sizes are 
respectively 0.0025 AU, 0.5°, and 0.002 for the eccentricity. The 
color scale is the stability index (logDi), and the level curves 

give the values (see Fig. |5]l. The cross marks the best fit 
SI, and the horizontal line in each panel, the orbits common to 
both maps (ie. the intersection of each map). The best fit lies 
very close to a stable resonant island where we can pick a stable 
solution with a only marginally higher than that of the best 
fit. As an example we picked such a solution (S2), marked by a 
white filled circle, by slightly increasing Oc from 2.4832 AU to 
2.49 AU. The complete set of orbital elements for S2 is given in 
Tab|2] 

3g2 = 0. As a consequence, one of them is superfluous. A new 
fundamental frequency associated to the resonance Ig replaces it. 

The solution S2 is trapped in the 5:1 mean motion resonance 
with the following main resonant argument: 

9 - Al,- 5Ac + nib + 3otc . (8) 

The variations of 6 versus time are plotted in Fig. [T] (green line) 
and exhibit a libration around 60 = 0°. We observe nonetheless 
that this libration is the modulation of several different terms 
with similar amplitudes of approximately 40°, but on different 
time scales. This leads to a libration with an amplitude that can 
be higher than 180°. 
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Table 3. Fundamental frequencies for S2. nt, and are the mean 
motions, gi and g2 are the secular frequencies associated with 
respectively TUt and Wc, and Ig is the libration frequency of the 
resonant angle 6 - Ab - 5 Ac + mb + 'i-ujc- 
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Fig. 7. Time variation of the resonant argument - Ah - 5 Ac + 
vjb + 3ro-c for the orbital solution S2 (green line). 6 is in libra- 
tion around 6*0 = 0°, with a modulation of two terms of period 
Pg - 80.13 years, and Pacj = 749.10 years, amplitudes of about 
35 and 50 degrees respectively. The black line shows the first 
secular term contribution. 



Table 4. Quasi-periodic decomposition of the resonant angle 
9 - Ab - 5Ac + TUb + 'i'UJc for an integration of the orbital so- 
lution S2 over 1 million years. The decomposition is given in 
the form given by Eq.|9] and each frequency Vj is expressed as a 
combination of the fundamental frequencies in the form given by 
Eq.[T3] For the sake of brievity, we only give the first 30 terms. 
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In order to describe more accurately the behavior of 0, we 
search for a quasi-periodic decomposition of 0{t). We st art with 
a fre quency decomposition using frequency analysis jLaskaH 
l2003i) as: 

0{t) = Y,AjCos(vjt + 4>j) (9) 

j 

And then we decompose each frequency vj on the four funda- 
mental frequencies {ric, g\, g2, k)' 

Vj = ajHc +I3jgi +6jg2 + Jjk , (10) 

where aj,/3j, 6j, and yj are integers. Each term of the decomposi- 
tion follows a D' Alembert-like relationship expressed in Eq.fTTl 

aj+^j + 6j = (11) 

which can be used to simplify the expression of vj. Indeed rear- 
ranging the right hand side of Eq.[TOl we can write: 

Vj = ajric +/3j(gi - g2) + (fij + 5j)g2 + Jjh , (12) 

and using Eq.[TT] we get: 

Vj = aj{nc-g2) +/3j{gi -g2) + 7jle ■ (13) 

The first thirty terms of the decomposition can be found in 
Table |4] The decrease of the amplitudes Aj is slow, due to the 
strong perturbations. 



The first term (j — 1) is responsible for the long term oscilla- 
tions with frequency -g2- The period coiTesponding to gi -g2, 
associated to the angle Am = Wb - Wc, is Pacj ~ 749.09 years. 
It is worth mentioning that the angle Atir is not in libration in 
this system, as opposed to what has been observed in planetary 
systems locked in a 1/2 mean motion resonance, such as GJ876 
(Laughlin & Chambers 2001; Lee & Peale 2002; Ji et al. 2001 
Cor reia et al.li2 009a) or HD89830 (iGozdziewski & Macieiewskil 
200 IT or in a 3/2 m ean motion resonance such as HD45364 
(ICorreia et al.ll2009bh . The second one (j - 2) introduces the 
short term oscillations of frequency Ig, corresponding to a period 
Pg ~ 101.49 years. More precisely the libration of is made of 
three difi'erent kinds of contributions: 

- secular terms: aj - jj - 0, hence of the form jSjigi — g2); 

- resonant terms: aj - and jj + 0, of the form - g2) + 

- and short period terms: Uj + 0. 

We plotted with a solid line in Fig.|7]the contribution of the term 
j - 2, the first secular contribution. 

4.4. Test particles 

We test in this section the possibilities for a third body around 
HD202206. To that end we integrate the S2 solution with added 
massless particles over 16000 years. As explained in the previ- 
ous sections, we compute two determinations n and n' of each 
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particles mean motion over two consecutive time intervals of 
8000 years. We then obtain a stability index D - log|« - n'\ 
for each particle, which we plotted with a color code in Fig. |8] 
Both panels are {a, e) grids of initial conditions of the test parti- 
cles. We vary the semi-major axis from 0.05 AU to 0.5 AU with 
a 0.0025 step size in the left panel, and from 0.5 AU to 10 AU 
with a 0.05 step size in the right panel. We span eccentricities 
of the test particles from to 0.9 with a 0.005 step size. Since 
the particles mean motions are higher in the left panel, they were 
integrated with a time step of lO^^* years instead of 2 x lO"-'. 

Due to the very large eccentricities of HD202206b and 
HD202206 c, the dynamical environment between and around 
them is very unstable. As a result, we don't expect any viable 
planets with a semi-major axis between approximately 0.12 AU 
and 6.5 AU (red dots). Most of these particles were actually 
lost before the end of the integration, either through collision or 
by having their eccentricity increased higher than 1 . The same 
computation with particles of one Earth mass yields very sim- 
ilar results. Assuming that S2 is a good representation of the 
HD202206 planetary system, we can use these results to put con- 
straints on hypothetical and yet undetected additional compan- 
ions. There is clearly two possible regions for new planets: either 
close to the star {a < 0.12 AU), or outside HD202206c {a > 6.5 
AU). 

In the first case, any planet massive enough should already 
have been detected as many full period are available in the data. 
Assuming a low eccentricity for the hypothetical companion and 
a 6 m/s instrumental precision we find that planets bigger than 
24 earth masses should have already been detected. A Neptune- 
sized planet can exists anywhere between 0.06 AU and 0.12 AU, 
and a 10 earth masses planet, anywhere between 0.02 AU and 
0.12 AU. 

In the second case {a > 6.5 AU) the period is greater than 16 
years, meaning that we have only covered approximately half an 
orbit at best. However a planet massive enough would create at 
least a detectable trend in the data. At 6.5 AU we can rule out 
the existence of a planet with more than half a Jupiter mass. At 
10 AU we can rule out a planet between 1 Mjup and 3 My„y„ de- 
pending on the phase. We conclude that a yet undetected planet 
smaller than half a Jupiter mass could exist at semi-major axis 
greater than 6.5 AU. 



4.5. Finding tine center of iibration inside a resonance 
ttirougli frequency anaiysis 

4.5.1 . Center of Iibration 

We will consider in this section the planar three-body problem, 
and more particularly, the planetary problem with a p : p + I 
mean-motion resonance. The problem is to find the orbit center 
of Iibration starting from a quasi-periodic orbit in the resonance. 

In the restricted case this center of Iibration is a well defined 
periodic orbit, such as the Lagrangian points in the 1 : 1 reso- 
nance. However in the general problem it is not so easy to define 
and to find this orbit. It has 4 degrees of freedom, and its quasi- 
periodic orbits live on 4-torus in the phase space. It should be 
noted that the problem can be restricted to 3 degrees of freedom, 
using the angular momentum reduction. Each dimension of the 
4-torus is associated to one of the four fundamental frequencies 
(/o>/i>/2,/3) of the orbit. Let's suppose that /o is the frequency 
of the resonant mode. The orbit equivalent to the center of Iibra- 
tion is living on a 3-torus, depending only on (/i , /2, /s). In other 
words, it is a torus where the fourth dimension associated to the 
resonance has a null amplitude. Intuitively, we can represent it 



Fig. 9. Schematic representation of the 3-torus of the center of 
Iibration. We represent the 4-torus of a given resonant orbit by a 
2-torus (a doughnut) as it is not possible to represent it otherwise. 
The center of Iibration 3-torus is then represented by the circle 
in the center of the interior of the doughnut. 



has the center of the 4-torus in the dimension associated with the 
resonance. 

4.5.2. Quasi-periodic decomposition 

If this orbit was periodic we could use a simple Newton algo- 
rithm to find it (provided that we start within the convergence 
radius), as it is a fixed point of a Poincare map. This method has 
been extensively used in numerical search of periodic orbit fam- 
ilies of the three body problem (see for instance Henon 1974], 
Il997h . 

We present here a new numerical method to find quasi- 
periodic center of Iibration, using the fact that we can get an ac- 
curate quasi-periodic decomposition of a numerically integrated 
quasi-periodic orbit with frequency map analysis ( Laskar 2003]) . 
Let X(f) be a state vector of such an orbit. Using frequency anal- 
ysis, we can obtain a quasi-periodic representation for any com- 
ponent X of the vector X: 



xit) = ^A(,)E'<^-^>' 

(k) 



(14) 



for each coordinate of the vector X, with k = (ko, k\ , k2, kj) e Z"*, 
/ - ifo, f\, fi, fi), and A(i) e C. We can separate these sums 
in two parts: x(t) = u{t) + v(t) where u does not depend on the 
resonant frequency /o, and v has all the terms depending on /q. 



'«(f)= 2 A(,,E'<*'/>' , 
v(/)= X A(,,E'<^-/>' . 



(15) 



The quasi-periodic orbit described by U(f) is precisely living 
on a 3-torus which has the characteristics we are looking for 
However it is probably not a solution of the equations of motion. 
But we can assume that it is close to one. Hence, we can use 
U(0) as a new initial condition, and obtain a new resonant quasi- 
periodic orbit X'(f) with: 



X'(0) = U(0) . 



(16) 



The amplitude of the resonant terms in this new orbit will be 
smaller. In other words, it lives on 4-torus closer to the 3-torus 
of the center of Iibration. We can then iterate this procedure to 
suppress all the terms with /o. 
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Fig. 8. Mean motion diffusion of tests particles, integration S2 solution with massless particles over 16000yr. We compute two 
determination n and «' of each particules mean motion over two consecutive time intervals of 8000 yr The color scale code the 
stability index log |« - n'\. Initial conditions for the massless particles are (/p, Mp, Wp, Qp) = (90°, 0°, 0°, 0°). On both panel we vary 
the initial eccentricity of the test particules from to 0.8 with a step size of 0.004. On the left panel the semi-major axis varies from 
0.05 AU to 0.5 AU with a step size of 0.0025 AU, and on the right panel it varies from 0.5 AU to 10 AU with a step size of 0.05. 
The white crosses mark the position of the two planets, and the collision lines are traced with white Unes. 



4.5.3. Application to HD202206 

We present here its application to the HD202206 system, using 
our orbital solution S2 as a starting resonant orbit. We work with 
the X - {ab,Ab,Zb,ac,^c,Zc) state vector, where Zh - e^E"^'' 
and Zc - CcE"^'. Each step p of the method is decomposed as 
follows: 

- Numerical integration of X'^''', 

- Determination of (/o, f\,f2,fi) using frequency analysis, 

- Quasi-periodic decomposition: x*'''(f) ^ u^P\t) + v*'''(f), 

- New initial conditions: x<''+"(0) = xW(0) - v('''(0) 

The initial conditions are not computed exactly following Eq.[T6l 
Since we work with a finite number of terms, and the amplitudes 
of the terms in v(f) are supposed to become small, the error we 
make because we only take into account that a finite number of 
terms will be smaller 

The convergence proved to be fast as we reduced the ampli- 
tude of the most important resonant terms by 2 orders of magni- 
tude in just 4 steps (Fig.fTTI). We show graphically the decreas- 
ing amplitude of libration at each step in Fig. [10] where we plot 
approximated sections of the successive torus projected in the 
iflc, TUc) plane. No resonant terms are found in the first 100 terms 
of the quasi-periodic decomposition of each variable at the last 
step (with the exception of Aj,). In fact, in half of the variables, 
there are no resonant terms left in the first 300 terms (Fig.fTTTi. 

Once we find an orbit with a zero amplitude libration mode 
to a good approximation, we also get the 3-Torus it is living on 
since its quasi-periodic decomposition gives us a parametriza- 
tion of the torus. The three angular variables of the torus appears 
in the decomposition: 

>1 = fit + ifri 
• (f>2 = fit + (1^2 (17) 
03 = + 



105 




2.5 2.504 2.508 2.512 2.516 2.52 2.524 
ao (AU) 

Fig. 10. Projection of a section of the 4-torus of each step's tra- 
jectory in the (uc, (jJc) plane. Each torus is approximated using a 
(truncated) quasi-periodic decomposition of the trajectory. 

where 1/^2, and ij/T, are initial phases. Let - E"^'. We can 
rewrite Eq.[T4]to reveal the underlying torus: 

401,02,03) = 2 ^wC^'^^r (18) 

W 

Of all the orbits living on this torus we can choose the closest 
to the radial velocity data. This can be done by minimizing the 
on the three initial phases 0i, tf/2, and 1^3. 

We note S3 the orbit obtained that way at step 6. We give 
initial conditions for S3 in table |5] This solution yields a square 
root of reduced equal to 1.55. If the system is locked in the 
5:1 mean motion resonance, the libration mode is likely to be 
dampened through dissipative processes. In that regard, the so- 
lution S2 is unlikely as it is on the edge of the resonant island, 
and exhibit high resonant mode amplitude (40° for the resonant 
critical angle). We expect that the real solution will be closer to 
S3. 
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Fig. 11. Evolution of the amplitude of the libration mode in the 
quasi-periodic decomposition of each variable at each step. In 
the top panel we plot the relative amplitude of the first term de- 
pending on /o compared to the first not constant term. In the 
bottom panel we plot the position of this first term in the decom- 
position. The first step (abscissa 0) is the S2 orbital solution, and 
the last step (abscissa 6) is the orbital solution S3 (table|5l). 

Table 5. Orbital parameters of an orbit close to the center of li- 
bration of the 5:1 mean motion resonance. This orbit was obtain 
from S2 (Table|2| at the sixth iteration. We expect that the actual 
system will be closer to this orbit as dissipative processes will 
have dampened the resonant mode. 
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5. Coplanar orbits 

5.1. Stability and low orbits 

In this section we investigate the system behavior when both 
planets remain in the same orbital plane (/ = 0°) and are pro- 
grade, but the inclination to the line of sight is lower than 90°: 

- ib = ic = i < 90° ; 

- Q, = 0° . 

For each inclination value, we compute a new best fit with a 
Levenberg-Marquard algorithm. The value obtained at dif- 



ferent inclinations is plotted in Fig. [T3] (solid line). Interestingly 
enough, we obtain better fits for lower inclination down to 15°. 

In this configuration, the masses of the two companions grow 
approximately in proportion with 1 / sin / when the inclination / 
diminishes. There is thus very little changes between 90° and 
50°, as seen in Fig. [12] (two top leftmost panels). For lower in- 
clinations, due to the increased masses, mutual perturbations be- 
come stronger, and less orbits are stable. However, amongst the 

lowest orbits, the ones in the 5:1 mean motion resonance 
remain stable for inclinations up to 15°. 

For / - 10°, no stable orbits are left (bottom rightmost 
panel). This puts clear limits on the inclination of the system, 
and on the masses of the two companions : m/, < 95.5 Mjup and 
mc < 12.5 Mjup. 

An interesting property showing in Fig. [12] is that the lowest 

orbits are always in the vicinity of the stable island of the 5:1 
mean motion resonance, for all inclinations. We believe that it is 
a strong point supporting the hypothesis of a HD202206 system 
locked in this resonance. It appears however, that after / = 50°, 
and for lower inclinations, the low region is slightly shifted 
towards lower values than the resonance island. 

In order to have a more precise picture, we span / values from 
10° to 90° with a step size of 5°. For each inclination value, we 
start by computing a new best fit with the Levenberg-Marquard 

algorithm. ™d the stability indexes Di and D2 are com- 
puted in the (oc, e^) plane of initial conditions, around the best 
fit. The size of each grid is 101 x 161 dots, with step sizes of 
respectively 0.0025 AU and 0.002. For each of these maps, we 
compute Dii,„ (see section [374] ), and detect stable regions, and 

their relative positions to the observations in terms of -^Jx^. 
In order to get a synthetic vision of all this data, we can 

get an estimation of the lowest stable ^/x^ for each inclination 

(Fig.[T3]l. We plot with respect to / the lowest ^/x^ of D, < Dn, 



^Xr 

orbits (broken curve), along with the best fit value (solid curve). 




Fig. 13. Evolution of the best fit y/xi (solid curve) and an esti- 
mation of the best stable orbits y/x^ in function of the inclination 
/. For each value of / between 90° and 10°, a new fit, assuming 
a coplanar configuration, is computed. We then look for stable 
orbits around each of these orbital solutions in a {Oc, e^) plane of 
initial conditions. 

It mainly confirms that the distance between stable orbits, 
and low -\lx^ orbits is really small from 90° to 50°. Is is actually 
slightly better between 60° and 50°. 
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Fig. 12. Dynamics of a coplanar HD202206 system for different value of the inclination /. Each panel is a diffusion map in the 
(Qc, Cc) plane of initial conditions constructed the same way as the top panel of Fig.|5](the top-left panel is in fact a copy). For each 
value of inclination a new fit, assuming a coplanar system, is computed. 




Fig. 14. Percentage of stable orbits inside a given 
curve against / (see Fig. [13]). 



level 



'X^ or- 



To have a better idea of the correlation between low 
bits and stable orbits, we look at the percentage of stable orbits 
inside a given ^jx^ level curves (Fig.fT4li. It gives a synthetic rep- 



resentation of the overlap between low ^]xi regions and stable 
regions. It is very clear that inclinations lower than 30°, although 
possible, are very unlikely. It also appears that inclinations be- 
tween 90° and 50° are the most probable. 

If we limit the acceptable inclinations to 30° < / < 90°, we 
can derive limits for the masses of the inner and outer planets 
(assuming a coplanar prograde configuration): 

- 1 < niblmf < 2 

- 1 < mjmf^ < 2 

It is interesting to estimate the time when the radial velocity 
data will allow to determine a more accurate estimation of the 
inclination /, and masses for the system. Assuming this model 
is close enough to the true system, we look at the differences of 
the radial velocities in the / - 90° case, and / = 50°, or / - 30°. 
With the hypothesis of an instrumental precision of 7 m/s, we 
find (Fig.fTSll that we will have to wait until about 2015 to sepa- 
rate between / - 90° and / - 30°, which corresponds to approx- 
imately a factor 2 in the masses. For / - 50°, five more years are 
needed for the differences between the two configurations to be 
greater than the measurements precision. 
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Fig. 15. Differences between the radial velocity of a stable so- 
lution with / = 90°, and / = 50° (red curve), or / = 30° 
(green curve). The starting epoch of the three integrations is 
2453000.00yZ). The blue horizontal lines represents the preci- 
sion of the cuiTent set of data, obtained with the CORALIE in- 
strument. The vertical blue line marks the beginning of 2009. 



5.2. Resonant and secular behavior 

We end this study of the coplanar configurations with a quick 
look at the dependence of the resonant and secular dynamics 
on the inclination. For each inclination, we picked a stable orbit 
with a low value, and plotted its fundamental frequencies 
le,gi,and §2 in Fig. [16] 
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Fig. 16. Libration and secular frequencies. For each inclination 
i we pick a stable orbit with a low and integrate it over 1 
million years. Using frequency analysis, we determine for each 
orbit its fundamental frequencies. The solid curves represent, for 
each frequency, a 1 / sin^ / law. 



As expected from the perturbation theory, when the masses 
increase, the secular frequencies also increase (in absolute 
value). We verify that it follows a rule in 1 / sin^ / (solid curve). 
This is a consequence of the fact that the most important terms 
responsible for the secular dynamics are of order two of the 
masses. 



6. Mutually inclined orbits 

In this section we drop the coplanarity constraint. We allow the 
inclinations ii, and ic to vary independently, and we also allow 
variations of Q.c, the longitude of ascending node of c. To span 
the possible values for ii,, ic, and Q^. in an efficient manner, we 
restrict ourselves to two mass ratios k/, (Eq.[3]) for planet b: 

- kh^l (ih = 90°), 

- kh^2 (ib = 30°), 

and three mass ratios k^ for planet c: 

- kc^l (i, = 90°), 

- kc^2 {ic = 30° or = 150°), 
_ kc^5{ic^ ll°or/c = 169°). 

For each couple of inclinations it, and ic, we let AO vary between 
0° and 360° with a 10° step size, and for each triplet {ih, ic, Af2) 
we perform a fit with the Levenberg-Marquard minimization, 
and compute a diffusion grid in the (Oc, ec) plane of initial con- 
ditions around the fit. The step sizes of the grids are respectively 
0.0025 AU and 0.004. We plotted in Fig. [HJ; and Fig. [HJ; the 
mutual inclination / as a function of AO for reference. For each 
configuration we look at the proportion of stable orbits inside 

the ^[)^ -1.5 level curve (FigfTTb andfTSb). We will assume 
that stable zones with ^[x^ < 1 .5 orbits harbor potential solu- 
tions of the system. Note that we keep the initial inclination //, 
in ]0; 7r/2] as opposed to ic (see section lTTT l. Also for any given 
values a e]0; 7r/2] and p e [0; 2k{, assuming all the other orbital 
elements identical, the configuration {ii, - njl, ic - a, D.c = /3) 
is symmetric to the configuration (ih = n/2, ic - 7T-a,Q.c = 

with respect to the plane ( i , k) (see Fig. |2]i. Since this plane 
contains the line of sight, the two configurations are indistin- 
guishable using radial velocity measurements. Hence for kb - I 
(it = 90°) we will not look at ic = 150° and 169°. 



6.1. kb = 1 

For ki, - 1 (ih = 90°, Fig. [TTT i we find that significant stable 
zones inside = 1-5 exist mostly for aligned and anti-aligned 
ascending nodes (i.e. AO ^ 0° and 180° in FigfTTb). 

For kc - 1 (red curves) the aligned configurations (AQ = 
0° + 10°) are coplanar and prograde (/ ^ 0°), and the anti-aligned 
configuration (AO - 180° + 10°) are coplanar and retrograde. 
Outside those two particular cases, we find no significant sta- 
ble zones for ^/x^ < 1.5. Indeed, while the resonant island is 
roughly centered on the lowest level curve, it is not stable 
outside the coplanar configurations. We find that there exists an 
extended zone of stability outside the resonance for AO ^ 30°, 
but it lies just outside the 1.5 level curve. Note that due to sym- 
metries, the situation for AO > 180° mirrors that of AO < 180°. 

For kc = 2 (green curves) we find again stable zones for 
AO ^ 0°, but not around 180°. However stable regions with 
potential solutions exist for AO values up to 90°, coiTesponding 
to mutual inclinations between 60° and 90°. This is mostly due 
to the minimum x^ getting smaller up to AO - 90° (see the 
green curve in Fig. [TTh). While the stable regions, both in the 
resonance and outside, shrink when AO augments, the x^ = ^-5 
level curves encompass a larger area. 

Finally for kc - 5 (blue curves), no significant stable zones 
are found at low x^ values. 

To summarize, potential solutions (stable orbit with a low 
X^) for ih = 90° mainly exists for coplanar configurations, inside 
the 5:1 mean motion resonance. If planet c's mass is kept low 
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Fig. 17. Stable configurations for ib - 90°, that is for the min- 
imum mass of planet b. For each value of ^ in (90°, 30°, 11°), 
and AQ between 0° and 360° we compute diffusion index over 
grids in the (qc, ee) plane of initial conditions with step sizes of 
respectively 0.0025 AU and 0.004. Each grid is also centered 
on the minimum computed for each couple (/c, AQ). We plot 
against AQ. the square root of the minimum in the top panel, 
and the proportion of stable orbits inside ^/x^ = 1 .5 in the mid- 
dle panel. The mutual inclination corresponding to each triplet 
(ib, ic, AQ) is plotted in the bottom panel. 



(kb < 5) stable regions do exist for non coplanar configurations, 
but they are located outside the lowest region. A noteworthy 
exception is the retrograde configuration, where stable orbit are 
only found close to the coplanar case, which only happens for 
kr = 1 and AQ ^ 180°. 



6.2. kb = 2 

When we double the mass of planet b (ib - 30°, FigfTsTl. once 
again retrograde potential solutions can only occur for copla- 
nar orbits: i^ = 150°, and AQ ^ 180° (green dotted curve). 
Concerning prograde orbits, there exists potential stable solu- 
tions for a mutual inclination of / ^ 60°: 

- i, = 90° and AQ = 0° (red curve), 

- ic - 30° and AQ = 180° (green solid curve). 



6.3. Conclusions 

To sum up we find that the configurations which have a signifi- 
cant stable zone at low^^ values are found mostly when the two 
lines of nodes are aligned. That is for AQ close to 0° or 180°. In 
addition, stable orbits with the lowest are all in the 5:1 mean 
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Fig. 18. Stable configurations for ib = 30°, that is for approx- 
imately twice the minimum mass: nib - 33.5 Mjup. See Fig. [TtI 
caption for more details. 



motion resonance except for nearly coplanar retrograde config- 
urations, where they can also be close to the commensurability. 
However we can also find stable orbits outside the resonance 



in the prograde resonance, usually with -\/xi higher than 1.45. 
Retrograde configurations seem to be limited to nearly copla- 
nar orbits with anti-aligned ascending nodes. Other than that, 
we could not find any clear correlation with mutual inclination. 

The fact that we find non-resonant stabl e solutions for retro - 
grade configurations is consistent with Smi th & Lissaueij(l2009l) . 
They have shown that retrograde configurations allow more 
closely packe d systems than prograd e configurations. It was also 
suggested by iGayon & BoisI (1200 8') that retrograde configura- 
tions are likely alternatives both from the radial velocity data 
and the long term stability point of view. Forming such a sys- 
tem remains however difficult, and we will thus not favor this 
hypothesis. 



7. Discussion and conclusion 

Assuming that the system is coplanar, we performed a system- 
atic study of the dynamics of the system for different inclinations 
to the line of sight. We are able to find constraints for the incli- 
nation to the line of sight: 30° < / < 90°. This means that the 
companions' masses are most likely not greater than twice their 
minimum values: 



< 2 



1 ^ / (0) 

- 1 < mblmy 

- 1 < mjmf^ < 2 



We also studied the influence of mutual inclination for two dif- 
ferent inclinations of the planet b (ib - 90° and //, = 30°), but did 
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not find any clear correlation other than that retrograde potential 
stable solution consistent with the radial velocity data seem to 
be limited to mutual inclinations close to 180° (i.e. nearly copla- 
nar orbits). As lGozdziewski et al.l(l2006l) . we find possible stable 
solutions with low;^'^ for a wide range of mutual inclinations be- 
tween 0° and 90°. The current data cannot yield more precise 
constraints. Also the masses determination is dependent on the 
stellar mass which is not well established. 

Although all published dynamical studies of HD202206 sug- 
gest that b and c are in a 5:1 m ean motion resonance, it is still a 
debated question. For instance. iLibert & HenrardI (l2007l) assume 
that it is just close to the commensurability. Libration of Act oc- 
curs for particular initial values of this angle, providing a sta- 
bilizing mechanism outside the mean-motion resonance, not far 
from the best fits. In most cases other than retrograde coplanar 
configurations, those orbits in near-commensurability are worse 
solutions than the ones in resonance, but they could be more 
probable if the eccentricities are overestimated (especially for 
b). We find that all significant stable zones with the best O-C are 
in the 5:1 mean motion resonance. In fact the minimum is al- 
most always in the resonance or very close to it, and stable orbits 
in the resonance can be found with not significantly higher 
than the best fit. In addition the O-C level curves tend to follow 
roughly the resonant island, even t hough the agreement is not as 
perfect as for the HD45364 system (Correi a et al.l2009b[) . This is 
an improvement from Correia et al. (2005) where the best fit lay 
outside the resonant island, and the;^^^ had to be degraded to find 
a stable solution. We thus believe that the resonant configuration 
is the most probable. We provide a stable solution (S2, Table|2]l 
in the coplanar edge-on case. This solution shows a high am- 
plitude resonant mode in the libration of the critical angle. We 
believe that this resonant mode is probably dampened by dis- 
sipative processes. We use frequency analysis to find a tore on 
which such orbits exist. Although the specific orbit we give in 

Table|5]does not have a very low at 1.55, we expect that the 
true orbit will be close to it with a low libration amplitude. 

Note that for retrograde configurations, the picture is quite 
different. The best fit lies in a very stable region just outside the 
mean motion resonance. While these orbits are valid candidates 
from the dynamical and the observational points of view, we do 
not favor them as the formation of these systems is hard to ex- 
plain. 

We investigated the possibility of undetected companions. 
We found that planets with masses smaller than approximately 
one Neptune mass can exist for semi-major axis lower than 0.12 
AU. 0.5 Mjup planets are also possible beyond 6.5 AU. No plan- 
ets are possible between 0.12 AU and 6.5 AU as they would be 
unstable. The two planets model may prove to be wrong in the 
future, but these hypothetical new companions should not have 
a big impact on the already detected ones. 
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